Defensive Resistance of Cowpea Vigna unguiculata Control Megalurothrips usitatus Mediated by Jasmonic Acid or Insect Damage

Vigna unguiculata is a vital vegetable crop in Southeast Asia, and Megalurothrips usitatus can cause huge damage to this crop. Enhancing the resistance of V. unguiculata against M. usitatus is a promising way to protect this crop; however, there is limited information regarding the mechanism underlying the resistance of V. unguiculata against M. usitatus. Here, a behavior assay was performed to explore the resistance of V. unguiculata against M. usitatus after insect damage or treatment by jasmonic acid (JA). Furthermore, transcriptome and metabonomics analysis was used to detect the putative mechanism underlying the resistance of V. unguiculata against M. usitatus. The pre-treatment of Vigna unguiculata with JA or infestation with Megalurothrips usitatus alleviated the damage resulting from the pest insect. We further identified differentially expressed genes and different metabolites involved in flavonoid biosynthesis and alpha-linolenic acid metabolism. Genes of chalcone reductase and shikimate O-hydroxycinnamoyltransferase involved in flavonoid biosynthesis, as well as lipoxygenase and acyl-CoA oxidase involved in alpha-linolenic acid metabolism, were upregulated in plants after herbivory or JA supplementation. The upregulation of these genes contributed to the high accumulation of metabolites involved in flavonoid biosynthesis and the alpha-linolenic acid metabolism pathway. These transcriptional and metabolite changes are potentially responsible for plant defense and a putative regulatory model is thus proposed to illustrate the cowpea defense mechanism against insect attack. Our study provides candidate targets for the breeding of varieties with resistance to insect herbivory by molecular technology.


Introduction
Megalurothrips usitatus (Bagrall) is a major pest worldwide and can reduce yield, causing significant economic losses each year [1]. M. usitatus usually destroys petals and fruits of plants, leading to severe yield failure [2]. Chemical pesticides, such as spinetoram, acetamiprid, chlorfenapyr, and fluridine nitrile, are frequently applied to control this pest; however, this rapidly causes severe resistance and environmental pollution [3]. Moreover, chemical control is usually difficult because this pest is able to hide inside the flowers and pods of plants [4]. Therefore, developing alternative approaches to control this pest is an urgent matter. Cowpea (Vigna unguiculata spp. sesquipedialis), which has high nutritional Table 1. Summary of insect feeding performance on plants. The statistical approach was the LSD test. The herbivory insect is M. usitatus, and the plant is V. unguita; 20 insects were placed on the surface of the leaf in each treatment group. The feeding index means the number of pest gathered on each leaf group/total number of insects.

Profiling of Transcriptome Sequencing and Gene Expression
To further explore the mechanisms underlying the insect resistance of plants induced by jasmonic acid and insect damage, transcriptome analysis of the leaves was performed. Because no high-quality reference genome was available, a full-length transcriptome was constructed. A total of 15,198,519 subreads were obtained with an average length of 2339 bp and an N50 of 2900 bp. After filtering out the low-quality reads (quality < 0.9 or full passes < 3), 275,887 full-length nonchimeric reads and 33,304 polished high-quality isoforms remained, and most isoforms ranged from 1000 to 5000 bp ( Figure 1A). Finally, 27,167 nonredundant transcripts were obtained after removing redundant reads. By map- ping these full-length transcripts to the NR, NT, Pfam, GO, and KOG databases, a total of 26,938 transcripts were annotated in at least one of these databases ( Figure 1B and Table S1). Healthy leaves (CK), leaves of herbivores after 24 h (24H) and 48 h (48H), and treatment with JA (JA) in three biological replicates were used for RNA-seq analysis. Ultimately, a total of 272,618,481 reads was generated. The total mapped read rates ranged from 73.84% to 76.68% (Table S2). Correlation analysis revealed that the correlations between the three replicates were reliable (Figure 2A). In the principal component analysis (PCA), the first principal component (PC1) explained 24.16% of the variation, and the second principal component (PC2) explained 13.67% of the variation. Moreover, samples were separated by JA treatment or herbivory attack along PC1 ( Figure 2B).
To further explore the mechanisms underlying the insect resistance of plants induced by jasmonic acid and insect damage, transcriptome analysis of the leaves was performed. Because no high-quality reference genome was available, a full-length transcriptome was constructed. A total of 15,198,519 subreads were obtained with an average length of 2339 bp and an N50 of 2900 bp. After filtering out the low-quality reads (quality <0.9 or full passes <3), 275,887 full-length nonchimeric reads and 33,304 polished high-quality isoforms remained, and most isoforms ranged from 1000 to 5000 bp ( Figure 1A). Finally, 27,167 nonredundant transcripts were obtained after removing redundant reads. By mapping these full-length transcripts to the NR, NT, Pfam, GO, and KOG databases, a total of 26,938 transcripts were annotated in at least one of these databases ( Figure 1B and Table  S1). Healthy leaves (CK), leaves of herbivores after 24 h (24H) and 48 h (48H), and treatment with JA (JA) in three biological replicates were used for RNA-seq analysis. Ultimately, a total of 272,618,481 reads was generated. The total mapped read rates ranged from 73.84% to 76.68% (Table S2). Correlation analysis revealed that the correlations between the three replicates were reliable (Figure 2A). In the principal component analysis (PCA), the first principal component (PC1) explained 24.16% of the variation, and the second principal component (PC2) explained 13.67% of the variation. Moreover, samples were separated by JA treatment or herbivory attack along PC1 ( Figure 2B).

Functional Classification of Differentially Expressed Genes Based on Kyoto Encyclopedia of Genes and Genomes Analyses
To dissect gene expression changes in plants after herbivory or being triggered by JA, differentially expressed genes (DEGs) were identified with the criteria of false discovery rate (FDR) <0.05 and absolute value of log2ratio ≥ 1. In the comparison of CK and 24H,

Functional Classification of Differentially Expressed Genes Based on Kyoto Encyclopedia of Genes and Genomes Analyses
To dissect gene expression changes in plants after herbivory or being triggered by JA, differentially expressed genes (DEGs) were identified with the criteria of false discovery rate (FDR) < 0.05 and absolute value of log2ratio ≥ 1. In the comparison of CK and 24H, 4043 DEGs (Table S3) were found, with 2311 downregulated and 1732 upregulated, and the number of DEGs increased to 5053 (Table S4) (2852 downregulated and 2201 upregulated) in CK vs. 48H; however, only 2567 DEGs (Table S5) (1327 downregulated and 1240 upregulated) were found in CK vs. JA ( Figure 3A). The number of upregulated DEGs was less than the number of downregulated DEGs in all the comparisons. Moreover, the number of DEGs in the herbivory group was higher than that in the JA treatment group, implying that gene expression levels changed more dramatically after insect attack than after JA treatment. A Venn diagram of DEGs was constructed to illustrate the three comparisons. A total of 1506 DEGs appeared simultaneously in all comparisons, and 1189, 1912, and 420 were unique to CK vs. 24H, CK vs. 48H, and CK vs. JA, respectively ( Figure 3B and Table S6).  Then, DEGs identified in each comparison were subjected to KEGG pathway analysis. Seventeen, sixteen, and fifteen pathways were found to be significantly enriched in CK vs. 24H (Table S7), CK vs. 48H (Table S8), and CK vs. JA (Table S9), respectively (Figure 4A-C). These pathways mainly consisted of amino acid, sugar, and secondary metabolism pathways. To more accurately confirm the common genes that occurred in all the comparison groups, KEGG analysis of common DEGs was conducted, and enriched pathways were related to phenylpropanoid biosynthesis, terpenoid backbone biosynthesis, flavonoid biosynthesis, alpha-linolenic acid metabolism, isoflavonoid biosynthesis, and biosynthesis of various plant secondary metabolites, among other pathways ( Figure 4D and Table S10). Then, DEGs identified in each comparison were subjected to KEGG pathway analysis. Seventeen, sixteen, and fifteen pathways were found to be significantly enriched in CK vs. 24H (Table S7), CK vs. 48H (Table S8), and CK vs. JA (Table S9), respectively ( Figure 4A-C). These pathways mainly consisted of amino acid, sugar, and secondary metabolism pathways. To more accurately confirm the common genes that occurred in all the comparison groups, KEGG analysis of common DEGs was conducted, and enriched pathways were related to phenylpropanoid biosynthesis, terpenoid backbone biosynthesis, flavonoid biosynthesis, alpha-linolenic acid metabolism, isoflavonoid biosynthesis, and biosynthesis of various plant secondary metabolites, among other pathways ( Figure 4D and Table S10).

Induction of Flavonoid Biosynthesis, Phenylpropanoid Biosynthesis and Alpha-Linolenic Acid Metabolism by Insect Feeding or Exogenous JA
Genes involved in flavonoid biosynthesis and alpha-linolenic acid metabolism have been proven to be responsible for resistance induced in plants (Table 2). Therefore, we further focused on DEGs associated with these pathways. Among the 24 DEGs involved in the phenylpropanoid biosynthesis pathway, most DEGs were upregulated, with seventeen and seven genes being up-and downregulated, respectively. The downregulated genes were two peroxidases [EC: 1.11

Induction of Flavonoid Biosynthesis, Phenylpropanoid Biosynthesis and Alpha-Linolenic Acid Metabolism by Insect Feeding or Exogenous JA
Genes involved in flavonoid biosynthesis and alpha-linolenic acid metabolism have been proven to be responsible for resistance induced in plants (Table 2). Therefore, we further focused on DEGs associated with these pathways. Among the 24 DEGs involved in the phenylpropanoid biosynthesis pathway, most DEGs were upregulated, with seventeen and seven genes being up-and downregulated, respectively. The downregulated genes were two peroxidases [EC: 1.11

Dynamic Transcriptome Responses to Insect Attack
From the K-MEAN analysis results, nine clusters representing 7717 genes are shown in Figure 5 (Table S11). Subclass 6 contained the most transcripts (2229), and subclass 9 contained only 406 transcripts. A total of 600 transcription factors belonging to 72 families were identified. In subclusters 2, 4, and 5, we found that the expression patterns in the treatment group were higher than those in the control group. However, clusters 3, 6, and 9 showed adverse trends. Among the TFs that occurred in subclusters 2, 3, 4, 5, 6, and 9, MYB (13), AP2/ERF-ERF (19), WRKY (14), and bZIP (12) were included. These TFs exhibit different treatment-specific expression patterns, and the specific expression patterns of these may reflect pivotal roles in plant defense.

Metabolite Profiles Indicate Differences in Metabolic Regulation under Herbivory Attack and JA Treatment
To further explore metabolite changes in plants, metabolomic profiles were generated for the control and treated plants. A total of 1224 metabolites were detected, and they were divided into 10 classes, among which flavonoids, phenolic acids, and lipids were main classes. Moreover, alkaloids, amino acids and derivatives, lignans and coumarins, nucleotides and derivatives, organic acids, phenolic acids, and terpenoids were also detected. The correlation analysis indicated good reproducibility of the different biologically replicated samples ( Figure 6B). In PCA, PC1 and PC2 explained 24.02% and 16.3% of the variation, respectively ( Figure 6A). PC1 separated the samples by herbivory attack or JA treatment. By analyzing the differentially expressed metabolites (DEMs), we found 259, 310, and 220 different metabolites in CK vs. 24H, CK vs. 48H, and CK vs. JA, respectively (Table 3 and Tables S12-S14). In the 24H group vs. CK, 24 metabolites were downregulated, and 205 were upregulated. In the comparison of the 48H group vs. CK, 66 and 244

Metabolite Profiles Indicate Differences in Metabolic Regulation under Herbivory Attack and JA Treatment
To further explore metabolite changes in plants, metabolomic profiles were generated for the control and treated plants. A total of 1224 metabolites were detected, and they were divided into 10 classes, among which flavonoids, phenolic acids, and lipids were main classes. Moreover, alkaloids, amino acids and derivatives, lignans and coumarins, nucleotides and derivatives, organic acids, phenolic acids, and terpenoids were also detected. The correlation analysis indicated good reproducibility of the different biologically replicated samples ( Figure 6B). In PCA, PC1 and PC2 explained 24.02% and 16.3% of the variation, respectively ( Figure 6A). PC1 separated the samples by herbivory attack or JA treatment. By analyzing the differentially expressed metabolites (DEMs), we found 259, 310, and 220 different metabolites in CK vs. 24H, CK vs. 48H, and CK vs. JA, respectively (Tables 3 and S12-S14). In the 24H group vs. CK, 24 metabolites were downregulated, and 205 were upregulated. In the comparison of the 48H group vs. CK, 66 and 244 were downor upregulated, respectively. Moreover, in the comparison of the JA and CK groups, 62 were downregulated and 158 were upregulated. The number of DEMs in CK vs. JA was less than those in CK vs. 24H and CK vs. 48H, which were similar to the number of DEGs in the transcriptome analysis. A total of 103 different metabolites were found to overlap in all three groups. The KEGG analysis showed that these different metabolites could be attributed to secondary and primary metabolite-related pathways, including flavonoid biosynthesis, alpha-linolenic acid metabolism, phenylpropanoid biosynthesis, and alanine, aspartate, and glutamate metabolism.   For better understanding of the changes in flavonoid and alpha-linolenic acid-related metabolites after herbivory or JA induction (Table 4), we identified the DEMs associated with flavonoid biosynthesis, alpha-linolenic acid metabolism, and phenylpropanoid biosynthesis, and found 13 DEMs. Metabolites associated with flavonoid biosynthesis included trans-5-O-(p-coumaroyl) shikimate, 3,5,7-trihydroxyflavanone (pinobanksin), dihydrokaempferol, naringenin chalcone, 5-O-p-coumaroylquinic acid*, liquiritigenin, naringenin, and isoliquiritigenin. DEMs involved in alpha-linolenic acid metabolism included 9-hydroperoxy-10E,12,15Z-octadecatrienoic acid, 9-hydroxy-10,12,15-octadecatrienoic acid, 12-oxo-phytodienoic acid, 13S-hydroxy-9Z,11E,15Z-octadecatrienoic acid, and jasmonic acid. All of these metabolites were induced after herbivory or JA treatment.   For better understanding of the changes in flavonoid and alpha-linolenic acid-related metabolites after herbivory or JA induction (Table 4), we identified the DEMs associated with flavonoid biosynthesis, alpha-linolenic acid metabolism, and phenylpropanoid biosynthesis, and found 13 DEMs. Metabolites associated with flavonoid biosynthesis included trans-5-O-(p-coumaroyl) shikimate, 3,5,7-trihydroxyflavanone (pinobanksin), dihydrokaempferol, naringenin chalcone, 5-O-p-coumaroylquinic acid*, liquiritigenin, naringenin, and isoliquiritigenin. DEMs involved in alpha-linolenic acid metabolism included 9-hydroperoxy-10E,12,15Z-octadecatrienoic acid, 9-hydroxy-10,12,15-octadecatrienoic acid, 12-oxo-phytodienoic acid, 13S-hydroxy-9Z,11E,15Z-octadecatrienoic acid, and jasmonic acid. All of these metabolites were induced after herbivory or JA treatment.

Discussion
JA and insect feeding were reported to play vital roles in plant-induced defense against herbivores because they can reduce insect damage [11,12]. This may be caused by the accumulation of defensive compounds induced by insect feeding or JA application. In our study, less damage was observed in plants pretreated by insect herbivory or JA than in healthy leaves. Drastic changes in the transcriptome and metabolome also occurred. In the group comparisons, 4043, 5053, and 2567 DEGs were identified in CK vs. 24H, CK vs. 48H, and CK vs. JA, respectively, and the metabolome analysis showed that there were 259, 310, and 220 different metabolites in CK vs. 24H, CK vs. 48H, and CK vs. JA, respectively. Generally, the insect-induced resistance of plant depends on producing a series of defensive compounds to affect the feeding, growth, and survival of herbivores [7][8][9]. In addition, the defense response of plants induced by JA or insect attack is mediated by these compounds in the phenylpropanoid biosynthesis pathway, flavonoid biosynthesis, and alpha-linolenic acid metabolism, which have been identified in many plants [13,14]. In our study, seven peroxidases [EC:1.11.1.7] involved in the phenylpropanoid biosynthesis pathway were found to be upregulated, which might play an important role in mediating the defensive responses of cowpea to insect damage [17] 1.195] involved in this pathway were also induced by insect attack or JA; these genes are known to perform critical functions under biotic stress condition [18]. In addition, two chalcone isomerases [EC:5.5.1.6], three chalcone reductases, and one shikimate O-hydroxycinnamoyltransferase [EC:2.3.1.133] were upregulated. These genes were associated with flavonoid biosynthesis. Flavonoids are the secondary compounds that play roles in plant defenses, which might be the result of changes in some metabolites that lead to the alleviation of insect damage. A similar study indicated that JA alleviates pest-associated losses in groundnut by mediating POD, PPO, and other defensive compounds, such as phenols, H 2 O 2 , and MDA (malondialdehyde) [19]. Plants protect themselves from pest attack with a large number of secondary metabolites [20]. These metabolites usually inhibit insect damage because they are feeding deterrents [21]. Plant metabolites may manipulate feeding behavior of herbivores [22]. Metabolites associated with flavonoid biosynthesis included trans-5-O-(p-coumaroyl) shikimate, 3,5,7-trihydroxyflavanone (pinobanksin), dihydrokaempferol, naringenin chalcone, 5-O-p-coumaroylquinic acid, liquiritigenin, naringenin, and isoliquiritigenin, which were herbivory-or JA-induced and accumulated in the cowpeas. Naringenin plays key roles in inhibiting the growth and development of the common cutworm Spodoptera litura [23]. The growth rate of tea geometrids showed significant decreases by feeding an artificial diet supplemented with naringenin [24]. In this study, seven metabolites and six transcripts in flavonoid biosynthesis were significantly correlated. The metabolites included trans-5-O-(p-coumaroyl) shikimate, pinobanksin, dihydrokaempferol, naringenin chalcone, liquiritigenin, naringenin, and isoliquiritigenin. The transcripts included chalcone isomerase [EC:5.5.1.6] (12h1_transcript_15683; 12h1_transcript_31775), chalcone reductase (12 h1_transcript_29690; 12h1_transcript_30450), and shikimate Ohydroxycinnamoyltransferase [EC:2.3. 1.133] (12h1_transcript_27132; 12h1_transcript_26689). Therefore, we assumed that DEGs would lead to a difference in secondary metabolite accumulation in cowpea, and then modulate the feeding behavior of thrips. A previous study confirmed that secondary metabolites are critical in modulating feeding behavior of insects [25]. The defensive roles of jasmonic acid and flavonoid pathways have been widely reported in a large range of plants [22].
Transcription factors are responsible for many biological processes in plants, including development, growth, and defense [34][35][36]. In our study, K-means clustering revealed that MYB (13), AP2/ERF-ERF (19), WRKY (14), and bZIP (12) showed different expression patterns in the pretreated plants, and CmMYB19, WRKY7, WRKY58, WRKY62, WRKY64, and WRKY76 were found to be induced under pathogen or insect attack [37,38]. The induction of bZIP TFs was observed to contribute to the defense response against greenbug infestation in sorghum [39]. Therefore, these transcription factors have a potential role in mediating plant defense.
In summary, the present study revealed that the pretreatment of Vigna unguiculata with JA or infestation with Megalurothrips usitatus contributed to induced insect resistance of host plants. The mechanism of insect resistance was associated with changes in secondary metabolites or phytohormone precursors at the transcriptome and metabolome levels, including chalcone reductase, chalcone isomerase, and lipoxygenase. Moreover, the defense mechanism involved a number of TFs such as MYB, WRKY, and bZIP ( Figure 8). Using molecular technology, the study provided candidate targets for the breeding of varieties with resistance to insect herbivory.

Plant Material, Cultivation, and Experimental Design of Megalurothrips usitatus Herbivory Treatment
The colony of M. usitatus used in this study was reared in the laboratory as described previously [3]. The cowpea V. unguiculata was planted in an artificial climate box with the relative humidity being 75%. Insects were reared on the cultivated cowpea in a growth chamber (26 ± 1 °C, 65 ± 5% relative humidity, 16:8 light/dark photoperiod). Plant materials were divided into three groups. The first group was infected with insects which were placed on the surface of plant leaves for 8, 16, 24, 48, and 72 h, respectively (including 4 treatments). The second group was sprayed with JA (including 1 treatment) and healthy plant leaves were used as the control group (including 1 treatment). The above three groups were then attacked by thrips. There were 6 treatments, namely, healthy control (CK), and 8, 16, 24, 48 and 72 h in this study. The leaves of these 6 treatments were used for the herbivory preference, and then the tested insect was used to choose the leaves from the 6 treatments. This process persisted for 48 h, and after that, the number of insects gathered on each specific treatment group was recorded, and the insect used in the choice experiment was also documented. The numbers of pests gathered on each leaf group/total number of insects were calculated using a stereoscope. A total of 20 M. usitatus for the insect damage was used in the experiment with three biological repetitions. Before being used, all the insects were not supplied with food for 24 h, and adult insects were used in this study.

Plant Material, Cultivation, and Experimental Design of Megalurothrips usitatus Herbivory Treatment
The colony of M. usitatus used in this study was reared in the laboratory as described previously [3]. The cowpea V. unguiculata was planted in an artificial climate box with the relative humidity being 75%. Insects were reared on the cultivated cowpea in a growth chamber (26 ± 1 • C, 65 ± 5% relative humidity, 16:8 light/dark photoperiod). Plant materials were divided into three groups. The first group was infected with insects which were placed on the surface of plant leaves for 8, 16, 24, 48, and 72 h, respectively (including 4 treatments). The second group was sprayed with JA (including 1 treatment) and healthy plant leaves were used as the control group (including 1 treatment). The above three groups were then attacked by thrips. There were 6 treatments, namely, healthy control (CK), and 8, 16, 24, 48 and 72 h in this study. The leaves of these 6 treatments were used for the herbivory preference, and then the tested insect was used to choose the leaves from the 6 treatments. This process persisted for 48 h, and after that, the number of insects gathered on each specific treatment group was recorded, and the insect used in the choice experiment was also documented. The numbers of pests gathered on each leaf group/total number of insects were calculated using a stereoscope. A total of 20 M. usitatus for the insect damage was used in the experiment with three biological repetitions. Before being used, all the insects were not supplied with food for 24 h, and adult insects were used in this study.

Gene Expression Analysis Using RNA-seq
Total RNA was extracted and purified using the Plant RNA Kit (BioTeke, Beijing, China) according to the manufacturer's instructions. The Agilent 2100 bioanalyzer was applied to assess the RNA integrity. The mRNA was enriched by magnetic beads with Oligo (dT) and the first strand of cDNA was synthesized using the fragmented mRNA as the template. Then, the RNA strand was degraded by RNaseH and the second strand of cDNA was synthesized using the DNA Polymerase I system. Subsequently, the cDNA library purity, quantity, and integrity were assessed using a Qubit 2.0 Fluorometer and Agilent 2100 Bioanalyzer prior to sequencing by an Illumina NovaSeq 6000 (Illumina Inc., San Diego, CA, USA).
Iso-Seq library construction was performed according to the Isoform Sequencing protocol (Iso-Seq) and the Clontech SMARTer PCR cDNA Synthesis Kit protocol (Mountain View, CA, USA). The sequencing data were processed using the SMRTlink software (Menlo Park, CA, USA). Briefly, subread BAM files were converted to circular consensus sequences (CCSs). Full-length nonchimeric (FLNC) sequences were selected from the CCS and then were subjected to nonredundant and cluster treatment by the ICE Quiver algorithm and Arrow polishing with nFL sequence. Finally, high-quality and polished full-length consensus sequences were produced [40,41].
Diamond blastx was used to perform gene function annotation with an E-value < 10 −5 [42]. The annotation datasets were taken from the NCBI nonredundant (NR) protein sequences, Swiss-Prot, Trembl, KEGG, GO, and KOG/COG databases. Transcripts were then aligned to the PFAM (Protein family) database to perform annotation using Hmmscan. Gene expression levels were calculated using RSEM [43]. FPKM was used to estimate gene expression levels. DEGs were analyzed using DESeq2 with cutoff padj < 0.05 and fold change > 2 [44]. GOseq and KOBAS were used for GO and KEGG enrichment analyses, respectively. The results of TF analysis of differentially expressed genes were extracted directly from the AnimalTFDB database [45].

Metabolite Profiling and Data Analysis
Biological samples were freeze-dried using a vacuum freeze-dryer (Scientz-100F, Anqing, China). Each freeze-dried sample was crushed using a mixer mill (MM 400, Retsch, Haan, Germany) with zirconia beads for 1.5 min at 30 Hz. Fifty milligrams of lyophilized powder were dissolved in 1.2 mL 70% methanol solution and vortexed for 30 s every 30 min for a total of 6 times. Following centrifugation at 12,000 rpm for 3 min, the extracts were filtered (SCAA-104, 0.22 µm pore size; ANPEL, Shanghai, China, http://www.anpel.com.cn/, access on 16 December 2022) before UPLC-MS/MS analysis.
Unsupervised PCA (principal component analysis) was performed by the statistics function prcomp within R (www.r-project.org, access on 16 December 2022). The Pearson correlation coefficients (PCCs) between samples were calculated by the cor function in R and presented as only heatmaps. For two-group analyses, differential metabolites were determined by VIP (VIP ≥ 1) and absolute Log2FC (|Log2FC| ≥ 1.0). Identified metabolites were annotated using the KEGG Compound database (http://www.kegg.jp/ kegg/compound/, access on 16 December 2022), and annotated metabolites were then mapped to the KEGG Pathway database (http://www.kegg.jp/kegg/pathway.html, access on 16 December 2022). Pathways with significantly regulated metabolites mapped to them were then fed into MSEA (metabolite set enrichment analysis), and their significance was determined by hypergeometric test p values.

Metabolomic and Transcriptomic Association Analysis
Pearson correlation coefficients were used to analyze the correlations between the differentially expressed genes from transcriptome analysis and the differential metabolites from metabolome analysis. If the correlation coefficient was less than 0, the correlation was negative, and vice versa. The correlated pairs with p values < 0.05 were used for further analysis. The differential genes and metabolites obtained were mapped to the KEGG pathway database simultaneously to obtain their common pathway information.